Load packages

library(LinQseq)

Step 1 Load data from 10X folder

# Load in the RNA UMI matrix
cbmc.rna <- as.sparse(x = read.csv(file = "../../Tutorial/data/GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz", 
    sep = ",", header = TRUE, row.names = 1))
cbmc.rna <- CollapseSpeciesExpressionMatrix(object = cbmc.rna)
# Load in the ADT UMI matrix
cbmc.adt <- as.sparse(x = read.csv(file = "../../Tutorial/data/GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz", 
    sep = ",", header = TRUE, row.names = 1))
cbmc.adt <- cbmc.adt[setdiff(x = rownames(x = cbmc.adt), y = c("CCR5", "CCR7", 
    "CD10")), ]

Step 2 Create object

# cbmc <- CreateSeuratObject(counts = cbmc.rna) cbmc[['ADT']] <-
# CreateAssayObject(counts = cbmc.adt)
cbmc <- createObject(rna = cbmc.rna, adt = cbmc.adt)

Step 3 Pre-process

1) Filter out unwanted cells (optional)

for this dataset, we don’t need to filter out unwanted cells

# cbmc <- subset(cbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 &
# percent.mt < softThreshold(cbmc))

2) Remove unwanted genes (optional)

for this dataset, we don’t need to filter out unwanted genes

# cbmc <- removeGene(object = cbmc,pattern = '^IG[HKL]')

3) Normalization

data Normalization for both ADT (CLR) and RNA (log)

cbmc <- dataNormalization(object = cbmc)

4) Indentify HVGs for RNA data

Call seurat function to identify highly variable genes (HVG) for RNA data

cbmc <- FindVariableFeatures(object = cbmc)  # directly use Seurat Function

5) Data scaling

Scale data for both ADT and RNA

cbmc <- dataScaling(object = cbmc)

Step 4 Linear dimension reduction (PCA)

directly call Seurat function for linear dimension reduction (PCA)

cbmc <- RunPCA(cbmc, features = VariableFeatures(object = cbmc), verbose = FALSE)  # directly use Seurat Function

Step 5 Determine number of PCs

call Seurat function JackStraw to determine number of PCs

# cbmc <- JackStraw(cbmc, num.replicate = 100) cbmc <- ScoreJackStraw(cbmc,
# dims = 1:20) JackStrawPlot(cbmc, dims = 1:15)
ElbowPlot(cbmc)

Step 6 Distance calculation and joint distance calculation

calculate cell-cell distances for RNA, ADT and joint. alpha was set to 0.5 as initial, number of PC was set to 20 by default.

cbmc <- jointDistance(object = cbmc)

Step 7 Non-linear dimension reduction (UMAP and t-SNE)

run UMAP as Non-linear dimension reduction for RNA, ADT and joint analysis.

# cbmc <- umapFromDistane(object = cbmc, assay = 'All')
cbmc <- tsneFromDistane(object = cbmc, assay = "All")

Step 8 Clustering

cbmc <- clusteringFromDistance(object = cbmc, assay = "All", resolution = c(0.2, 
    0.2, 0.2))

Step 9 Visualization ADT vs RNA vs Joint

1) Cell clusters

# gridDimPlot(cbmc, wide.rel = 1.5, legend = FALSE, reduction.prefix =
# 'tsne_', height.rel = 0.5)

plots <- generateGridDimPlot(cbmc, legend = FALSE, reduction.prefix = "tsne_", 
    darkTheme = FALSE)
listPlot(object = plots)

###### user also can only plot some of those plots by index, figure ident or
###### figure map info listPlot(object = plots, fig.ident = 'RNA')
###### listPlot(object = plots, fig.ident = 'RNA', fig.map = 'RNA')

###### user can use plotInfo() function to get index, figure ident and figure map
###### information, then plot figures by index
plotInfo(plots)
# listPlot(object = plots, fig.id = 1)

2) Heat maps

# Heatmap for joint clusters
heatMapPlot(object = cbmc, group.by = "jointClusterID", height.rel = 1, adt.label = TRUE)

# Heatmap for RNA clusters
heatMapPlot(object = cbmc, group.by = "rnaClusterID", height.rel = 1, adt.label = TRUE)

# Heatmap for ADT clusters
heatMapPlot(object = cbmc, group.by = "adtClusterID", height.rel = 1, adt.label = TRUE)

Step 10 Build trajectory using MST method

# Build MST
cbmc <- buildMST(object = cbmc, assay = "Joint")
# Build PG
cbmc <- buildPG(object = cbmc, assay = "Joint", reduction.prefix = "tsne_")
# Pseudotime
cbmc <- pseudoTime(object = cbmc, method = "MST", group.by = "jointClusterID", 
    graph = "jointTree", dist = "jointDist", root = "3")

Step 10 MST Trajectory Visualization

p1 <- trajectoryMSTPlot(object = cbmc, group.by = "jointClusterID", reduction = "tsne_joint", 
    pg = "jointPG")
p2 <- pseudoTimeMSTPlot(object = cbmc, reduction = "tsne_joint", pg = "jointPG", 
    pseudotime = "jointTreeMSTTime")
plot_grid(p1, p2, ncol = 2)

Step 11 build trajectory using kNN method

# Build KNN
cbmc <- buildKNN(object = cbmc, assay = "Joint", k = 10)
# Pseudotime
cbmc <- pseudoTime(object = cbmc, method = "KNN", group.by = "jointClusterID", 
    graph = "jointGraph", dist = "jointDist", root = 3)

Step 11kNN Trajectory Visualization

p1 <- trajectoryKNNPlot(object = cbmc, group.by = "jointClusterID", reduction = "tsne_joint", 
    graph = "jointGraph", line.size = 0.1)
p2 <- clusterConnectionPlot(object = cbmc, group.by = "jointClusterID", reduction = "tsne_joint", 
    graph = "jointGraph", line.color = "black", cutoff = 0)
p3 <- pseudoTimeKNNPlot(object = cbmc, reduction = "tsne_joint", graph = "jointGraph", 
    pseudotime = "jointGraphKNNTime", line.size = 0.1)
plot_grid(p1, p2, p3, ncol = 2)